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NUMERICAL APPROXIMATION OF LEVEL SET POWER MEAN 

CURVATURE FLOW 

AXEL KRONER*, EVA KRONER*, AND HEIKO KRONER* 


Abstract. In this paper we investigate the numerical approximation of a variant of the mean 
curvature flow. We consider the evolution of hypersurfaces with normal speed given by H k , k > 1, 
where H denotes the mean curvature. We use a level set formulation of this flow and discretize the 
regularized level set equation with finite elements. In a previous paper we proved an a priori estimate 
for the approximation error between the finite element solution and the solution of the original level 
set equation. We obtained an upper bound for this error which is polynomial in the discretization 
parameter and the reciprocal regularization parameter. The aim of the present paper is the numerical 
study of the behavior of the evolution and the numerical verification of certain convergence rates. 
We restrict the consideration to the case that the level set function depends on two variables, i.e. 
the moving hypersurfaces are curves. Furthermore, we confirm for specific initial curves and different 
values of k that the flow improves the isoperimetrical deficit. 

Key words, geometric evolution equations, level set formulation, viscosity solution, finite 
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1. Introduction. Geometric evolution equations, especially curvature-dependent 
interface motion has been studied for many years in both pure and applied mathe¬ 
matics. One application example is the evolution of soap films and the behavior of 
the boundaries of oil drops on a surface of water which evolve into disks. In mate¬ 
rial science, for example, the evolving surfaces might be grain boundaries in alloys 
which separate differing orientations of the same crystalline phase. In image pro¬ 
cessing, for example, one wants to identify a dark shape in a light background in a 
two-dimensional image. Therefore a so-called snake contour is evolved so that it wraps 
around the shape. We refer to [23j for a more detailed exposition of these applications 
and, e.g., [231 CDS HU HQ and references therein for further applications. 

The most famous case for such an interface motion is the mean curvature flow 
of closed n-dimensional hypersurfaces in Euclidean space R n+1 or as special case the 
curve shortening flow of closed curves in the plane. Under this flow the hypersurface 
moves in normal direction so that the normal speed equals the (mean) curvature and 
a convex initial hypersurface shrinks to a round point in finite time, where in case of 
closed, embedded plane curves this even holds if the convexity assumption is left out, 
cf. [37] and [29] ED [36] . 

Mean curvature flow can be formulated in parametric form, where the moving 
hypersurface is given by a parametrization over a fixed hypersurface which depends 
on the evolution time as variable, cf. a special case is graphical mean curvature 
flow, where the hypersurface is given as the graph of a height function, cf. [55]. A third 
possibility is a phase field approach to mean curvature flow, cf. [44] and references 
therein, and the fourth which will be considered in the following in more detail is to 
consider the PDE resulting from the level set formulation. The level set formulation is 
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powerful because it can handle topological changes of the moving hypersurface. Level 
set methods were introduced by Sethian and Osher, see gaEUES], and have been 
applied to a wide range of problems. 

In the present paper we are concerned with the level set formulation of a modified 
version of the mean curvature flow. Instead of the mean curvature we prescribe the 
normal speed of the evolution to be a power k > 1 of the mean curvature and assume 
that the initial hypersurface has positive mean curvature. In our previous paper gO] 
we proved a priori error estimates for a finite element approximation of this flow. The 
aim of the present paper is the numerical study of the behavior of the evolution and 
the numerical verification of certain convergence rates. We restrict the consideration 
to the case that the level set function depends on two variables, i.e. the moving 
hypersurfaces are curves. Furthermore, we confirm for specific initial curves and 
different values of k that the flow improves the ’isoperimetrical deficit’. 

To specify how our flow looks like we give the following parametric formulation 
of this flow. Let M be a smooth n-dimensional compact manifold without boundary 
(at the moment it is sufficient to assume only k > 0) and Xg : M —> R n+1 a smooth 
embedding such that Mq = xg (M) has positive mean curvature, then we consider a 
solution of the following fully nonlinear parabolic initial value problem. Find T > 0 
and a smooth mapping 


(1.1) 

x : [0, T) x M —► R n 

with 



*(0,-) = x 0 , 

(1.2) 

d 


s *(f,£) = -ffV 


Here, H and v denote the mean curvature and the outer normal of x(t, M) at x(t, £), 
respectively. We call this a power mean curvature flow (PMCF). 

Why is this flow interesting and how does this flow behave? This flow has been 
considered in a series of papers under different aspects. In [50] it is shown that the 
flow (1.2) exists on a maximal, finite time interval and that, approaching the final 
time, the surfaces contract to a point. In m the flow is considered in the case k > 1. 
It is shown that if initially the ratio of the biggest and smallest principal curvature 
at every point is close enough to 1, depending only on k and the dimension n of the 
hypersurfaces, then this is maintained under the flow. As a consequence the authors 
of [511 obtain that, when rescaling appropriately as the flow contracts to a point, the 
evolving surfaces converge to the unit sphere. The paper [45] shows that for k > n — 1 
the flow improves a certain ’isoperimetrical difference’. As singularities may develop 
before the volume goes to zero, a weak level-set formulation for such flows is developed 
and it is shown that the nronotonicity of the isoperimetrical difference is still valid. 
This proves the isoperimetrical inequality for n < 7. A further reason which makes 
this flow interesting is that mean curvature flow is used in denoising images, cf. [T] for 
such applications, and we expect that PMCF with the possibility to choose different 
values for k > 1 is an interesting alternative for this purpose. We mention the remarks 
in [43] which state that in the case 0 < k < 1 equation (1.2) plays a key role in the 
context of image processing. 

As announced above we will perform our calculation for curves, so it is of interest 
to give references for this case, both of theoretical and practical nature. In mm a 
finite element approximation of the parametric formulation of the flow (1.2) in the 
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case of curves is formulated and stability bounds are derived, see also [Sj. In [42] the 
evolution of plane curves driven by a nonlinear function of curvature and anisotropy 
is considered with a focus on the analysis of the parametric formulation of such a flow, 
see also EH- The analysis of boundaries of shapes in the context of morphological 
and shape image processing leads to an equation of the form (1.2) in the case of 


curves. This has been introduced in |H EU48] and we mention especially the case 


k = | which is the so-called affine curvature equation, cf. [5] and see also (1.8) and 


the following text lines. In summary one can say that our flow (1.2) (in the case 
k > 1) plays an important role for applications and has received a lot attention so far 
but the numerical approximation of its level set formulation has apart from our own 
previous work [40j not been analyzed yet. 

We introduce our notation. The Euclidean norm of R n is denoted by | • |. For 
an open subset f 1 of R” and m £ N*, p > 1 we denote the corresponding Sobolev 
spaces by W m ’ p (fl), W 0 m ’ p (fi), H m (D) = W m - 2 (fl) and i? 0 m (fl) = W 0 m ’ 2 (ft). The 
dual spaces are denoted by W~ m ' p {M) = IE™’ p (fI)* and the dual pairing by 


(1.3) 


w~ m ’ p {n) x w 0 m>p (fi) 9 (f, 1-9 (f, tp) = Ftps 


In the following we will introduce a (stationary) level set formulation for 0> 
and start for this purpose by recalling the (time-dependent) level set formulation for 
the mean curvature flow, i.e. the case k = 1 in (1.2). Let Mq C R n+1 be a given 
initial hypersurface and choose a continuous function u o : R n+1 —► R such that 


(1.4) 

If u : [0, oo) x R” +1 

(1.5) 


M 0 = {x& R n+1 : ii 0 (a:) = 0}. 

> R is the unique viscosity solution of 

jt u=]Du ' dW {m) =H{Du ’ D2u) 


in R n+1 x (0, oo) with u(0, •) = Uq in R" +1 , where 


( 1 . 6 ) 


h(p,y) = Y,( s ‘>- P ^:) Yu 


for p = (pi) € R" \ {0} and Y = (Yy) e R nxn . We call the family of the 
(1.7) M(t) = {x € R" +1 : u(t,x) =0}, t> 0, 


a (time dependent) level set mean curvature flow. Equation (1.5) is a quasilinear 


degenerate and possibly singular (if Du = 0) parabolic equation. Existence and 
uniqueness of a solution for this equation is proved in m hhi m- 

If we include our nonlinearity H k of the mean curvature in this formulation we 


get instead of (1.5) the fully nonlinear, degenerate and possibly singular parabolic 
equation 


(1.8) 4.-U*| (£(T-5^i D,D,U 


dt 


1,3 


\Di 




\Du\ 


1 -k 


To our knowledge an existence proof for (1.8) is only known for the case 0 < k < 1, 
cf. |43| and the references therein. But the case under consideration in the present 
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paper is k > 1. In case k > 1 the proof presented in [43j does not work any more 
because the linear growth of the elliptic part of the operator needed to apply classical 
arguments is not available. 

in the k = | case, i.e. the affine curvature 


The time dependent formulation 


1.8 


equation, is used for image processing, c'. [U(3Q]. In [43] equation ( |1.8[ ) in case 0 < 
k < 1 is approximated by a family of regularized equations and rates of convergence 
of the corresponding solutions are obtained. 

In the case k = 1, of course, we have mean curvature flow, and the corresponding 
equation ( |1,5 ) has been studied intensively analytically and numerically, cf., e.g., 
mmm 39 . We want to point out the paper m by Deckelnick, where the solution 
u e of a regularized version of (1.5) is approximated by a finite difference scheme which 
was originally proposed by Crandall and Lions [2DJ. In Deckelnick’s paper rates for 
the convergence of the discrete solution to the solution u of the (not regularized) level 
set equation are proved. The total error consists of a regularization error of the form 


(1.9) 


\u-U e \\ L oc in) < C a £° 


with a £ (0, 5 ) arbitrary and c a a positive constant, see [21], Theorem 1.2] for details, 
and a discretization error which is a polynomial expression in the numerical parameter 
and the reciprocal regularization parameter. Furthermore, the concrete value for the 
convergence order of the discretization error (and hence for the total approximation 
error) is very low; the main point here is that this rate is of polynomial order. 

This is not self-evident as can be seen in the paper (23- There the viscosity 


solution u of (1.51 is approximated by a solution u e of the regularized equation and 


then the regularized equation is approximated by a solution u Ei h of a semi discrete 
problem. The regularization error is again of the form ( |1.9[ ) but the error u e — u e ^h 
measured in a certain energy norm, cf. Theorem 6.4], is only of order c £ ft., where, 
and this is the important point, the constant c e depends exponentially on K Numer¬ 
ical tests as written there, however, suggest that the resulting bound overestimates 
the error. In the special case of two dimensions, i.e. the moving hypersurfaces are 
curves, Deckelnick and Dziuk JM] prove L°°-convergence (without rates) of the dis¬ 
crete solution provided h = h{e) sufficiently small, where ’sufficiently small’ is not 
given by an explicit formula or polynomial dependence. 

Let us now consider the case k > 1 which is the relevant one in the present paper. 
To circumvent the above mentioned problem with the growth in the elliptic part of 
the operator if k > 1 Schulze |5D[ uses a stationary level set formulation at which 
the nonlinearity due to the exponent k affects only lower order terms. We present 
Schulze’s stationary level set formulation of the PMCF. 

Let fl C R n+1 be open, connected and bounded having smooth boundary dfl 
with positive mean curvature. Here, dfl plays the role of the initial hypersurface. 
We call the level sets F t = d{x £ SI : u{x) > t}, t > 0, of the continuous function 
0 < u £ C°(S2) a (stationary) level set PMCF, if u is a viscosity solution of 


( 1 . 10 ) 


div ( 

\\Du\J 


1 


\Du\i 


u =0 


in fl 


on 9S7. 


For a definition of a viscosity solution for this equation we refer to [401 Section 2]. 
If u is smooth in a neighborhood of x £ S2 with non vanishing gradient and satisfies 


in this neighborhood (1.101, then the level set {u = u(x)\x £ S2} moves locally at 
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x according to (1.2). Using elliptic regularization of level set PMCF we obtain the 
equation 


( 1 . 11 ) 


div 


Du e 


E 1 + \Du E \ 2 


= -{£ + \Du e Y) 




U e = 0 


in Q, 
on dCl, 


which has a unique smooth solution u e for sufficiently small e > 0, cf. ESI Section 4]; 
moreover, there is cq > 0 such that 


( 1 . 12 ) 

uniformly in £ and (for a subsequence) 


C 1 ( n) < co 


(1.13) 


u e ^ u £ C 0 ’ 1 ^) 


in C°(f2). We call u a weak solution of (1.10), which is unique for n < 6. All the 
above facts are proved in E3 Section 4] under the assumption that k > 1. A weak 
solution of (1.10) satisfies (1.10) in the viscosity sense, cf. Section 113 Section 2]. 
Furthermore Schulze’s existence result is restricted to the case k > 1. 

What looks as a disadvantage at first glance, namely the fact that our level set 
function does not depend on the time explicitly and hence no explicit Euler method 
is applicable (as for example in Deckelnick’s paper mi), has the advantage that we 
have a divergence structure for the elliptic part of the operator which would not be 


the case if we would use the time dependent level set formulation (1.8) (in addition 
we would lack a proof of existence of a solution). 

To our knowledge our previous paper |40] is the only numerical analysis result 
for Schulze’s level set formulation so far. In [321 we proved an explicit rate for the 


convergence of the solution u e of (1.11) to the solution it of equation (1.10) which 


depends on fc, see Sect ion|3j where we recall the result. Using the divergence structure 
of the elliptic part of (1.11) we proved existence of a finite element approximation u e h 


of u e and an approximation rate. Summarized we have a total approximation error 


(1.14) 


u-u e h = (u-u e ) + ( u e - u e h ) 


which consists of the regularization error (first bracket in equation (1.14)) and the 


discretization error (second bracket in equation (1.14)). We obtained a polynomial rate 
in 1 and h for the total approximation error provided the discretization parameter is 
sufficiently small compared with the regularization parameter (the coupling between 
the discretization and the regularization parameter is also of polynomial order). The 
order of convergence is polynomial (in contrast, e.g., to [23[ Theorem 6.4], where the 
authors obtain exponential order of convergence) and comparable to the one obtained 
in m- In both cases despite from being of polynomial order the precise order is 
rather of theoretical value. As in the remarks following [551 Theorem 6.4] stating 
that experiments indicate that the proven approximation error overestimates the real 
error, we have a similar behavior for our situation, cf. Sections [111 and [4] 

Furthermore, we validate for some examples that the flow improves (i.e. decreases) 
the isoperimetrical deficit 


(1.15) 


A(f)"™ -Cn+lV(t), 
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where A(t) denotes the n-dimensional volume of the evolving hypersurface, V(t) the 
(n + l)-dimensional volume of the enclosed subset of R n+1 and c n+ 1 the Euclidean 
isoperimetrical constant, cf. [i50] for a proof of this property. The isoperimetrical 


deficit (1.15) is nonnegative and zero if and only if the hypersurface is a sphere. Since 


we restrict ourselves to the case n = 1 the evolving hypersurfaces are curves and the 
variables in the isoperimetrical deficit become arc length and enclosed area. And there 
holds C 2 = 47t. Geometrically more interesting is the case n > 2 because then one 
can have non convex initial hypersurfaces with positive mean curvature which develop 
topological changes under the flow. To validate our level set ansatz with respect to 
its convergence behavior and numerical properties restricting ourselves to the curve 
case seems to be a reasonable and legitimate technical simplification. 

The paper is organized as follows. In Section [2] we prove an error estimate for the 
discretization error and validate it numerically. In Section [3] we recall and calculate 
in detail an error estimate for the regularization error from our previous paper [40] 
and provide numerical examples. In Section [4] we provide numerical examples for the 
total approximation error. In Section [5] we illustrate with an example the influence of 
the exponent k on the behavior of the flow. 

2. Discretization error. In this section we assume that the space dimension 
n + 1 is 2 or 3 and that is convex. The latter is only a restriction if n + 1 = 3 
since dfl has positive mean curvature by our assumptions in the introduction. We 


fix e > 0 at a small value and approximate the solution u e of equation (1.11) by a 


finite element solution u e h which seems to be an appropriate method in view of the 
divergence structure of the operator. The goal is to analyze the discretization error 
u e - u e h . 

Let (24, lift,) be a quasi-uniform triangulation of Q with mesh size 0 < h < ho, ho 
sufficiently small, and 14, C H l (VLh) the finite element space given by 


( 2 . 1 ) 


V h = {u £ C°(n h ) : v\ dQh = 0, v\ T linear VT £ T h } . 


In view of the convexity of fl there holds f4 C fi. A function Uh £ 14 will be also 
considered as a function on fl by extending it by zero in fl \ Q/, . Then Vh £ H 1 ^). 
Our variational formulation is given as in our previous paper [40j by 


( 2 . 2 ) 


(Du e h , Dv h ) 
yje 1 + \Du E h \‘ 


dx = (e 2 + \Du e h \ 2 ) ™v h dx \/v h €V h . 

J 


For formal reason we might consider boundary tetrahedrons (boundary triangles in 
case d = 2) to be extended to a boundary tetrahedron with one ’curved face’. There¬ 
fore we will replace a boundary element T £ X4 (i.e. n + 1 vertexes of T lie on <9f 1) 
by T = T U B with 


(2.3) 


B = {tp + (1 — t)Pp | 0 < t < l,p £ F}, 


where F is the boundary face of T, i.e. n + 1 vertexes of F lie on dCl, and Pp is the 
unique minimizer of dist(p, -)| 9 n- We denote the resulting triangulation by Th■ This 
leaves the space of finite element functions we use (namely 14) unchanged. Note, that 
the boundary strip fl \ has measure 0{h 2 ). 


A similar equation as (2.2) is the stationary level set formulation for the inverse 
mean curvature flow which is used in [28j . There also a total approximation error, a 
discretization error and a regularization error appears. Furthermore, a rate for the 
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discretization error ( 0(h ) for the iJ 1 -error and 0(h 2 ) for the L 2 -error) is proved. But 
the dependence of the constants on the regularization parameter which appear in these 
error estimates is not analyzed theoretically. In contrast to m in gl3 no theoretical 
estimate for the regularization error (and hence for the total approximation error) is 
given, and such a rate seems to be an open problem so far, cf. |28j Remark 4]. But this 
issue is addressed numerically in [28] and calculations suggest that the regularization 
error is O(e), where £ here also denotes the regularization parameter. 

We remark that when we considered the discretization in our previous paper, 
see [301 Section 6], the space V),. consisted of continuous functions on ft^ which are 
piecewise polynomials of degree < 2 (and not linear as here) and we assumed that 
ft C R 2 and existence of a solution of (2.21 in this case was shown. The reason for this 


is that in EDI Section 6] error bounds for the discretization error are proved which 
contain the dependence of e explicitly. Therefore an estimate for the norm of the 
inverse of L e and its dual L *, where L e is the derivative of the regularized differential 
operator, see ( 2.81 for a definition, is calculated via the intermediate step of some 


rather technical sup-norm estimates and inverse inequalities which make it necessary 
to consider higher order elements. Since these £-dependencies do not play a role for 
the discretization error under consideration in this section we present the proof for 
the present case in easier form and also for W 1,p -norms in Theorem 2.2 with general 
p > n + 1 (contrary to Section 6], where p < 4 is assumed) which is necessary to 
prove an optimal L 2 -error estimate, cf. Theorem 2.3 


We start with a definition and properties of the linear operator L t and its dual. 
Let p > 1. We define for £ > 0 and z £ R” 


(2.4) 


:= fe(z) := 


and denote derivatives of f e with respect to z* by D z > f e . There holds 


(2.5) 


Zi „ „ „ / n $ij ZiZj 


D z if e {z) = |—f-, i>,iiW e (*) = rn --nf 


We define the operator 4> e by 


(2.6) 4> £ : W^(Q) W- 1 *' (ft), * e («) = -Di ( 

\\ Dv UJ \Dv\1 


where ^ = 1, so that (1.11) can be written as 


(2.7) 


$ e (u £ ) = 0. 


We denote the derivative of $£ in u e by 

(2.8) L e := D$ e (u e ) 

and have for all ip £ WQ ,p (ft) that 


(2.9) 


L e ip 


= ~D i (D z iD zj f e (Du e )D j ip) + ^-f e (Du e ) 1 * D zi f e (Du s )Djip 
=: -Di(a v Djip) + b'Diip, 


where we use the convention to sum over repeated indices. The coefficients a lJ and 
b l are in C 00 (ft). Note, that the estimate (1.121 is not available for higher order 
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derivatives of u e but since we fix e in the present section, this does not have an effect 
on the following considerations. 

The linear operator 


( 2 . 10 ) 


L e : W 0 1,p (f2) W- 1 - P *(n) 


and its adjoint operator L* are topological isomorphism, 
Appendix. 

We get from |T3] Theorem 8.5.3] for L = L e or L = L* 
there is a unique solution Uh G 14 of 


cf. Corollary 7.2 in the 


and F GW that 


( 2 . 11 ) 


(Lu h , ip h ) = Fip h \/ip h G 14, 


where u G H 1 (ft) is the unique solution of Lu = F and there holds the estimate 


(2.12) ||'U/i||iv 1 'P(fi) + ||u - u^||w 1 >p(n) < c||u||w 1 .p(n)- 

Furthermore, if F G L P (Q) we have 

(2.13) ||u — Uh\\w 1 ^(Q.) + h\\u — Uh\\Lp(n) < ch 2 \\F\\ LP {Q,y 


Remark 2.1. Note, that we actually used the assertion of [151. Theorem 8.5.3] 
with slightly different assumptions. The difference from the assumptions we need to 
the one assumed in m Theorem 8.5.3] are as follows. 

(i) We assume a right-hand side F G W~ 1,p (f 1) (instead F G L p { 12)). 

(ii) We consider the equation on (instead of a polygonal domain) and use as 

discretization the triple ( 24 , Q. 14). 

There holds the following theorem. 

Theorem 2.2. For every p > n + 1 and small h > 0 there exists a constant 
0 < c = c(||it e 


w 2 ’ 2 (Q)iP) such that (2.2) has a solution u £ h G 14 satisfying 


(2.14) 


- ujj||tyi,p(Q) < ch. 


This solution is unique in a small W 1,p -neighborhood of u e in 14. 

Proof. The proof of this lemma is adapted from [40] Section 6], where the result 
is proved for p < 4 and quadratic finite elements. 

We set 


(2-15) — {v h G 14 : ||u“ — Vh\\w 1 -p(n) < p}; 

where we choose 


(2.16) p=h x 

for an arbitrary and now fixed < A < 1. 

We will obtain Uh as the unique fixed point in B £ of the operator T : 14 —> 14 
with 


(2.17) L s (w h - Tw h ) = $ s (w h ), w h G 14- 

We show that B^ 0, that T is a contraction and that T(Bp) C B^. 
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(i) Let IhU £ be the interpolation of u £ , i.e. the continuous piecewise linear function 
on f \ which is equal to u £ at all nodes of fWe extend I^u 6 by zero to a function 
on Q. In view of 

(2.18) \\IhU e — w e ||u/ 1 ’P(n) < ch 

we have IhU £ G Bp for small h. 

(ii) Let v h ,w h £ Bp, £ h = v h - w h then using 


2.17) we conclude 


(2.19) 


L e (Tv h - Tw h ) = L e £ h + $ e (»h) - $ E (v h ) 

= {l s - D$ e (v h + e^))a 

=: F 


with a 0 G (0,1). In order to estimate \\F || w _i, P * (II) which leads to an estimate 
of ||Tvh — Twh\\w^<v(p) in view of Corollary 7.2 we choose £ Wq’ p (tt) with 
IIV’llw 1 ^* (n) < 1 and estimate (F,ip) . To do so we use a mean value theorem for 
which we need the following auxiliary estimate 


( 2 . 20 ) 


\Du £ -(Dv h + QDZ h )\\ La o {n) 

< || Du" — DI h u £ ||i=o(n) + || DI h u £ - Dvh\\L°°(Q) 

_ n+1 

< ch +cph p , 


where Vh = Vf,. + ©£h G Bp and where we used an inverse estimate. The resulting 
estimate implies 


( 2 . 21 ) 


\\Tv h -Tw h \\ w i, P tm <c(h +ph p 


W 1 ’ P(O) 


1 

<- 

“4 


w l -p(n ) 


for small h. 

(iii) Let Wh G B 1 }. There holds 


( 2 . 22 ) 


\\Twh—u £ \\ w ^,p(Q.) 

<|| Twh — TIhU £ \\w 1 -p(Q) + || TIhU e — LhM e ||vE 1 'P(n) 
+ || Ihu £ — M £ ||w 1 ’P(n) 

<^ + || TIhU E — /? i u e ||vE 1 .p(n) + ch 


It remains to estimate the norm on the right-hand side. There holds 

|| TIhU e — 7/iM E ||w 1 .p(n) <c||$ £ (/?iM e )|| w- 1 'P*( n) 

(2.23) =c\\<b e (Ii l u £ ) — 3 , e('a e )||w- 1 'P*(o) 

<ch 


again by a mean value theorem estimate. In view of (2.16) there holds 
(2.24) 

□ 


T(B h p ) C B h p . 
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In the following theorem we improve the L p -error estimate of Theorem |2. 2 [ there¬ 
fore we use a duality argument as in [28] . 

Theorem 2.3. For p > n+ 1 there holds 

(2.25) \\u E - u E h \\ LP{n) < ch 2 

with c = c(\\u e \\ W 2 , 2 ^,p) > 0. 

Proof. From the definitions of u E and u E h we get for all iph G 14 

(2 ' 26) In (iW|7 " ' WVhdX + In ( |V “‘ 1 *’ “ |V “ SI '‘) ' fhdx = ° 

This equation can be written as 


(2.27) 

with 


{AhVe e h )-Vip h dx+ [ (a E h -Ve E h )ip h dx = 0 

! JCl 


Al = 


= f D 2 f £ (\7u £ + t\7(u £ h — u £ ))dt 

Jo 

z h=\J fe(^u E + tV(u E h -u E ))*~ 1 Df e (S7u E + tS7(u E h -u E ))dt 


tth= k 

e h = u h - ^ 

and for later purposes we set 

A% =D 2 / £ (V« e ) 


(2.29) 




We define <p £ Wl' p (f2) by 
(2.30) L* e tp = \el\P~ 1 sgn{e e h ) 

and let tph £ Vh be the finite element solution of this equation. We test (2.301 with 
e e h and get in view of the symmetry of A e h that 

f K\ P dx = [ C A h^e s h ) • Vip h dx + f {a E h - Ve E h )ip h dx 

Jfl JQ 

[ {{A E h - A E h )Ve E h ) ■ Vip h dx + f ((4 - Oh)' V e%) ip h dx 
Jo. Jfl 

<c I |Ve|| 2 |V ^|dx + c I \Ve E h \ 2 tp h dx 

Jo. Jn 


(2.31) 


In view of Corollary 7.2 and (2.30) we get 


(2.32) 


\w 


l fc.|lw 1 >j>*(fi) Sc 


<c( / dxY = ellefIlf 


In 


c hllLp(n)> 
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so that 


(2.33) 


\e-h\\Lp(Q.) < ch 2 


Here, we used (2.311 and Theorem 12.21 to estimate \\e E h \\ w 1 ,^( 0 ,)- □ 


In the following we validate the convergence rates of Theorem 1 2.2 1 and Theorem |2.3| 
with a numerical example. Figure [l] shows the discretization error in the case of a 
unit circle as initial curve and e = 0.1 fixed. The discrete solutions for different 
values of the discretization parameter h are compared with the discrete solution u e 
on a fine grid with grid size h = 0.005. We compare on fio.oos and extend u e h = 0, 
where it is not defined. The discretization errors ||u e — ||n h for || • ||n h = || • IIl 2 ^)) 


= 


I H^Qh) 


and 


I 


Theorem 12.31 and Theorem 12.21 


|L°°(f 2 h ) are plotted and behave as shown in 



Fig. 1 : Discretization error for the unit circle as initial curve where k = 1, e = 0.1. 


In Figure[2]we see that the discretization error close to the boundary has a pattern 
which results from the approximation of the smooth domain Q by the polygonal 
domain The discrete solution u e h is zero at the polygonal boundary while u e 
is zero on the curved boundary dfl. 

3. Regularization error. We specify our a priori estimate for the regularization 
error presented in [2D] which depends on the choice of certain constants. Let 

(3.1) 7 > 1 + k 

and a, s > 0 so that 


(3.2) 
where 

(3.3) 


/3i(a, s) > /? 2 (a,s), 


Pi(a,s) 


2 — s + a(2 — \) 

7(2-i) + i-i’ 


02 (a,s) 


a + ks 
7 — k — 1 


0 < r < 


a 

7* 


and choose 
(3.4) 
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Numerical approximation of level set power mean curvature flow 



Fig. 2 : Discretization error for the unit circle as initial curve (we compare with a discrete 
solution on a fine grid with size ho). The pictures show the difference uf l0 — for e = 0.09, 
k = 1, ho = 0.025 and a) h = 0.4, b) h = 0.2, c) h = 0.1, d) h = 0.05. 


Note, that ( |3.2| ) obviously holds for sufficiently small a , s. There holds the following 
theorem, cf. [401 Theorem 3.1]. 

Theorem 3.1. There is c = c(k, fl) > 0 such that 


(3.5) 


K - w|| C 0(ft) < C£ n 


i(r,s) 


for all e > 0. 

The order of convergence stated in the previous theorem can be written more 
explicitly which is content of the following lemma. 

Corollary 3.2. For the evolution with normal speed H k , k > 1, the regulariza¬ 
tion error with respect to the C°-norm is of order O(ex) for all A > 2k. 

Proof. We rewrite the right-hand side of the estimate in Theorem 3.1 In (3.2 1 
we may assume that a, s are related by s = “. We multiply (3.2 1 by 7 and get 


(3.6) 


2 — s + a (2 — j:) a + sk 
2 


Now, we multiply with the denominator of the left-hand side and sort by a and s on 
each side which leads to 


(3.7) 


2 — s + a[2 — — > a 


1 _ k±l 

1 


*-i + *(2-i+ 7 
1 -rrm- 

7 


and after rearranging terms to 


(3.8) 


2 >s 


fc+i 

7 


+ k ( 2 ~ k + v) + ( 2 ~ + 1 ) 


1 _ fc+1 

7 


We may let 7 tend to infinity without changing the value of s (by adapting a corre¬ 
spondingly). Hence the right-hand side of (3.8) converges to 4/cs as 7 —> 00 and the 
claim follows. □ 
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We recall an interpolation lemma, cf |33; PDE II, Lemma 1.4.13]. 
Lemma 3.3. For 0 < /3 < a < 1 and a function v : 12 —> R holds 

(3-9) [v\p < 2 1 _ “Ms||v||co(q), 


where these expressions might become oo and 

IK®) ^ v{y) | 


(3.10) 


W " = S !*-»!■ 


Since u e is uniformly bounded in the (K-norm, cf. (1.12), we can use Lemma 3.3 to 
get 


(3.11) 


-U e \\ C 0,S ( Q)<c{/3)£ X{1 t5) 


for every 0 < /3 < 1 and 0 < A < ^. 

In the case k = 1 which means mean curvature flow we can realize in Corollary |3.2| 
every power of e which lies in (0, f). This is in accordance with the corresponding rate 
for the case of the time dependent level set regularization as considered in Deckelnick’s 
paper [2TJ Theorem 1.2] and Mitake’s paper [43] Theorem 1 ], 

The following numerical examples indicate that the regularization error is even 
smaller than stated in Corollary |3.2[ As a first example we calculate the regularization 
error in the case of the evolution of a unit circle as initial curve for which the exact 
solution of equation (1.10) is known. Let dB r 


,(0) C 


i.e. a circle with radius 


r 0 > 0 , be the initial curve then the exact solution u is given as 


(3.12) 


u(r) = 


W+l r k +1 
'0 ' 

k + 1 ’ 


where r denotes the radius variable in polar coordinates in M 2 with center in 0. As 
special case we choose k = 1 , r 0 = 1 , i.e. 

(3.13) u(r) = * - T —. 

In Figure [3] the error ||u e — u|| is plotted in this special case (and for k = 1.5 and 
k = 2), where || • || stands for || • || = |HU 2 (Q h ), INI = IHI/fi(W) and || • || = || ■ |U°o (fJh ). 
We remark that our theoretical estimate in Corollary |3. 2 1 does not provide information 
about an estimate with respect to || • ||jyi(n h )- The functions u e are calculated by using 
linear finite elements on a fine grid with mesh size h = 0.0125. The L 2 -error converges 
a little bit faster and the H^-e rror a little bit slower than of quadratic order to zero. 

In Figure [4] the scenario is the same as in Figure [3] apart from the fact that we 
now choose an ellipse (half axes lengths 1 and 2) as initial curve. Furthermore since 
we do not have an exact solution u for this case we use instead a solution u e h with 
e = 0.1 and small h = 0.0125. In Figure [5]we plot a section (along the long and short 
half axes of the initial curve) of the solution u e in the case of the circle and in Figure 
[ 6 ] in the case of an ellipse as initial curve for different values of e. 

In accordance with our a priori estimate in Corollary |3.2| the regularization error in 
Figure [3] seems to become larger for increasing k. This can be also seen from Figure[5] 
There we also observe that the approximation quality around the singularity of the 
flow deteriorates for k = 2 when changing £ from 0.5 to 0.17. 
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Numerical approximation of level set power mean curvature flow 





e 


Fig. 3: Regularization error in case of a circle as initial curve for k = 1,1.5, 2. 


Figure [7] shows level sets of u 01 for the case of the ellipse as initial curve and 
different values of k. We remark that our theory covers only the case k > 1 but in the 
special case of convex curves we also have a level set solution for general k > | which 
follows from the classification of the behavior of the evolution of curves by powers 
of the curvatures presented in Section 1 of [?j. Our observations are as follows. For 
k = 0.5 we see for £ = 0.1 a quite good approximation of the phenomenon of shrinking 
to a ’round point’ and further lessening of e does not show significant improvements. 
For all k the inner level line for e = 0.1 seems to be already ’round’ while for k = 2 
this seems to be far from a ’point’. 


4. Total approximation error. In the inequalities (2.14) and (2.25) appear 
constants c on the right-hand sides which depend on the solution u e of the regularized 
equation. To get an estimate for u — u e h in terms of e and h one has to make this 
dependence explicit. In our paper [40] we showed that there is a 7 > 0 such that 
if we couple h and £ by h = e 7 and use finite elements of order 2 (and quadratic 
boundary approximations), then there holds that the error u — u e h converges to zero 
with a polynomial rate in h with respect to the sup-norm. To estimate with respect 
to the sup-norm is natural since u is (only, in general,) a C -limit of u e and (as 
viscosity solution) Lipschitz continuous. A value for 7 and the convergence rate can 
be obtained by adapting it at each stage of the proofs in [40] as described there which 
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e e 



E 


Fig. 4: Regularization error in case of an ellipse as initial curve for k = 1,1.5, 2. 


leads to a rather technical large value of no practical interest. The main point is that 
we have a polynomial rate and not an exponential rate. As said before and explained 
by comparing the situation with I2TT Theorem 6.4], where the authors proved even 
only an exponential estimate, this is non-trivial. We let us furthermore inspire from 
the scenario of ]2J ;j Theorem 6.4] which overestimates the error rate as practical 
results indicate, cf. 23]. Therefore we start our calculations with the from practical 
point of view comfortable setting of continuous and piecewise linear finite elements, a 
polygonal boundary approximation and a coupling between e and h by setting e = h 
which already lead to convergence. 

Figure [8] shows the total approximation errors for the unit circle as initial curve in 
the cases k = 1,1.5, 2. Although we have only for the sup-norm a theoretical estimate 
we also plot the iJ 1 -error; we remark that in the situation of the circle the solution 
is of class C°°(B ro ( 0)). Figure [ 9 ] shows the same scenario as Figure [8] apart from the 
fact that we now consider the ellipse (half axes with lengths 1 and 2) as initial curve. 
Furthermore, as reference solution we consider a solution with h = 0.05 and e = 0.05. 

5. Effect of k on behavior of the flow for an example case. The phe¬ 
nomenon of becoming round can be measured by the isoperimetrical deficit 


(5.1) 


/(f) 2 — 47ra(t), 
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Numerical approximation of level set power mean curvature flow 



Fig. 5 : Radial solution for a circle as initial curve for k = 1, 1.5, 2. 


where l(t) denotes the length of the curve and a(t) the enclosed area at time t. 
According to theoretical results in E03 we confirm the monotonicity of this deficit 
during the evolution in the special case of the ellipse as initial curve, see Figure |10[ 
Furthermore, we see that with increasing k (and e = 0.05) the curves transform faster 
into a circle (they are not yet shrinked to a point except for k = 0.5, see Figure[7j where 
the ’point’ is reached quite well). In Figure[5]and Figure[6]we see when comparing the 
exact solutions for the circle for different values of k and the approximate solutions 
for the ellipse with e = 0.15 for different values of k , respectively, that the flow reaches 
the singularity for larger k earlier. 


6. Implementation. To calculate the finite element approximation of u E we 
used a discretization with unstructured grids, see Figure[lTJ These were generated by 
the mesh generator Gmsh, see [33]. We solved the non-linear equation (2.2) with a 
Newton method which uses a bi-conjugate gradient stabilized solver (BiCGSTAB) and 
SSOR preconditioning. For the implementation we used PDELab, a discretization 
module for solving PDEs which depends on the Distributed and Unified Numerics 
Environment (DUNE). As further references concerning PDELab we refer to [47],ITl]. 
information about DUNE can be found in [121H ES US]. In order to get solutions 
for small e we used a warm-start, i.e. we decreased e stepwise to the desired small 
value and performed on each stage a calculation with the solution for the previous e 
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0.1 


0.01 


0.001 L - 
0.01 



k = l.O k = 2.0 

1 04 L0*L2 


Fig. 7: Solution for ellipse for e = 0.1. 



Fig. 8: Total approximation error, e = h, for k — 1,1.5, 2, in case of a circle as initial curve. 
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0.01 0.1 1 

e 

Fig. 9: Total approximation error, e = h. for k = 1,1.5, 2 in case of ellipse with half axes 1 
and 2 as initial curve. 


as initial value. 

7. Appendix. Since L e : Hq(SY) — > iJ _1 (fl) is a topological isomorphism by 
classical L 2 -theory this also holds for L* : H ( \ (0) —> iJ _1 (12). We define the to L e 
associated uniformly, elliptic regular Dirichlet bilinear form of order 1 by 


(7.1) 
and set 

(7.2) 


B : WQ ,p (fl) x Wq' p (12) —> R, B[u,v] = f a^DiuDjV + VDiuv dx 

J n 

N p . ={w £ Wq’ p (12) : B[i/j,v] = 0 for every if) £ C'“(r2)} 

N p ={v £ W 0 1>p (12) : B[v,4>] = 0 for every </> £ 


From Fredholm’s alternative, cf. [54} Theorem 10.7], we deduce that for every F £ 
W _1,p (12) the equation 


(7.3) B[u, y\ = Fip M<p £ Wo’ p ’ (11) 

has a solution u £ Wq’ p (12 ) if and only if 


(7.4) 


v £ N p * => Fv = 0. 
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Numerical approximation of level set power mean curvature flow 



Fig. 10: Isoperimetrical deficit. 



Fig. 11: Mesh for the discretization with size h = 0.15 for ellipse with half axes 1 and 2. 


If dim N p * = dim N p = 0 then for every F £ W 1,p *(f2) equation (7.31 has a unique 


solution. 


Lemma 7.1. dim N p * = dim N p = 0. 


Proof. Let v £ N p *. From [54] Theorem 7.6] we get v £ W^ ' p ffl ) for all 1 < 
p 1 < oo, especially for p’ = 2. Since we know from L 2 -theory that (7.3) has a unique 
solution u £ Wg 1,2 (f]) if p = 2 and F = 0 we deduce that v = 0. Analogously we 
obtain the remaining claim. □ 

By bounded inverse theorem we conclude the following result. 

Corollary 7.2. L e ,L* are topological isomorphisms. 
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